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ABSTRACT 

This  paper  is  concerned  with  Lax-Wendroff  methods  for  a  class  of 
hyperbolic  history  value  problems.  These  problems  have  the  feature  that 
globally  (in  time)  smooth  solutions  exist  if  the  data  are  sufficiently  small 

'  fkt. 

and  that  solutions  develop  singularities  for  large  data.  prove  (second 
order)  convergence  of  the  Lax-Wendrof f  method  for  smooth  solutions  and 
investigate  numerically  the  dependence  on  the  initial  data,  we  demonstrate 
the  occurrence  of  shock  type  singularities  and  compare  the  results  to  the 
quasilinear  wave  equation  (without  Volterra  term). 
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SIGNIFICANCE  AND  EXPLANATION 

The  notion  of  viscoelastic  materials  can  be  modelled  by  partial 
integrrodif ferential  equations.  For  several  model  problems,  recent 

investigations  have  been  concerned  with  the  question  whether  or  not  these  ^  ^ 

equations  allow  the  development  of  shocks.;  In  this  paper,  we  investigate  this  w 
question  numerically.  Convergence  proofs  for  a  Lax-Wendroff  type  method  are 
given.  Our  computed  solutions  confirm  the  predictions  of  previous  papers. 

For  small  data,  the  solutions  remain  smooth,  but  for  large  data 
discontinuities  develop. 


The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive 
summary  lies  with  MRC,  and  not  with  the  authors  of  this  report. 


LAX-WEMDROFF  METHODS  FOR  HYPERBOLIC  HISTORY  VALUE  PROBLEMS 
Peter  Markowich*  and  Michaal  Renardy** 

1.  INTRODUCTION 

Viacoalaatic  materials  are  modelled  by  constitutive  lam  relating  the  stress  to  the 
history  of  the  strain  [6],  [14],  [15].  For  moat  of  the  constitutive  laws  suggested  by 
rheologists  (for  reviews  see  e.g.  [1],  [2),  [1C],  [20]),  the  functional  relationship  has 
the  fora  of  a  convolution  integral.  In  many  cases,  the  resulting  integrodifferential 
equation  of  aotion  can  be  regarded  as  a  perturbation  of  a  hyperbolic  equation  [17]. 

This  paper  deals  with  the  nuaerical  analysis  of  a  model  equation  of  this  fora  for  one- 
diaensional  viscoelastic  solids,  introduced  by  Daferaos  and  Nohel  [3] ,  [4] .  The  equation 
has  the  fora 


(1.1) 


tt 


t 

♦<u  )  -  /  b(t 

XX  Q 


s)f(ux(x,s)  )x<Js  +  f(x,t)  , 


where  b  is  a  positive,  bounded,  smooth,  integrable  kernel  and  are  smooth  functions 

satisfying 

«• 

♦(0)  -  *(0)  -  0,  ♦■(0)  >  0,  (0)  >  0,  ♦'(0)  -  **(0)  /  b(  T)dT  >  0  . 

0 

Particular  interest  in  this  problea  arises  for  the  following  reasons  It  is  well  known  that 
the  quasi linear  wave  equation 
(1.2)  »tt  -  ♦(«,)„ 

need  not  have  globally  smooth  solutions  even  if  the  initial  data  u(t  «■  0),  ufc(t  «  0)  are 
smooth.  Also,  no  dashing  occurs  since  the  energy 


(1.3) 


B[u] 


/  u^(x,t)dx  +  /  4(ux(x,t))dx 
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u 

is  a  constant  of  the  notion.  Hare  #(u)  -  /  4(o)do.  on  the  other  hand,  Dafemos  and 

0 

Hotael  13],  14)  have  shown  that  solutions  to  (1.1)  on  a  finite  interval  (with  Dirichlet  or 
Neunann  conditions  on  the  boundary)  remain  snooth  and  decay  to  sero  if  the  initial  data  and 
the  forcin9  tern  are  snooth  and  snail  in  appropriate  Sobolev  norns. 

The  dissipative  influence  of  the  integral  tern  thus  acts  against  the  formation  of 
singularities.  However,  this  effect  is  only  strong  enough  for  snail  initial  data.  The 
analysis  of  similar  equations  [8],  (13],  [24]  has  shown  that,  for  bounded  integral  kernels 
(this  is  essential,  cf .  [18] )  and  large  data,  smoothness  can  be  lost  in  finite  tine  and 
shocks  (i.e.  discontinuities  in  ux  and  ufc)  can  develop. 

In  this  paper  we  discuss  La*-Wendroff  type  schemes  for  the  numerical  solution  of 
(1.1).  The  equation  is,  for  this  purpose,  transformed  to  a  system  by  setting  ux  “  v, 
u^  *  w. 

In  section  3  we  show  that  -  for  appropriate  integration  methods  approximating  the 
Volterra  ten  -  these  schemes  converge  of  second  order  on  any  finite  tine  interval  on  which 
a  smooth  solution  exists.  For  the  stability  of  this  method,  the  nondegeneracy  condition 
♦'  *  0  is  essential,  as  is  evident  fron  recent  work  of  Friedel  and  osher  [5]. 

In  section  4  we  show  convergence  uniformly  up  to  t  •  •  for  the  case  of  spatially 
periodic  small  data  and  the  special  class  of  kernels 

M  -1,0 

(1.4)  b(0)  -  l  X  e  i  X,,  1,  >  0  . 

(■t  *  *  * 

Such  kernels  are  connonly  used  in  rheology.  The  proof  is  patterned  after  a  new  proof  for 
the  existence  of  globally  snooth  solutions,  which  we  present  in  section  2.  The  essential 
ingredients  for  this  proof  are  the  stability  of  the  trivial  solution  u  -  0  and  the 
quasllinear  hyperbolic  nature  of  the  equations.  If  the  1^  are  very  big,  a  stiffness 
problem  arises  in  the  nuaerical  analysis,  and  we  point  out  a  simple  modification  of  our 
scheme  which  avoids  this. 
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Section  5  is  concerned  with  numerical  experiments.  Computations  demonstrate  that 
singularities  occurring  for  large  data  are  indeed  of  ehock-type  (▼  and  w  have  jumps). 
Since  the  Imx-Wendroff  method  is  an  artificial  viscosity  method  the  shocks  are  not  sharp. 

A  "shock  layer"  occurs  whose  width  is  proportional  to  the  meah  size  (see  also  Kreiss  [12]). 

Our  computed  solutions  to  (1.1)  are  compared  to  those  for  the  corresponding 
quasilinear  wave  equation  (1.2).  This  comparison  shows  that  the  dissipative  mechanism  of 
the  Volterra  term  delays  the  time  of  the  blow-up  even  if  it  is  not  strong  enough  to  avoid 
the  developaMnt  of  singularities. 

It  is  well  known  that  -  for  a  scalar  conservation  law  -  the  Lax-Wendroff  method  may 
converge  to  a  nonphysical  weak  solution  (which  does  not  fulfill  the  entropy  condition),  see 
[5],  [7].  For  the  problem  (1.1),  no  theory  of  weak  solutions  and  entropy  conditions  has 
been  developed  yet,  and  the  main  interest  here  is  in  computing  globally  smooth  solutions  or 
smooth  solutions  up  to  the  blow-up  time.  For  this  purpose  the  second  order  Lax-Wendroff 
scheme  is  superior  to  first  order  methods,  which  do  converge  to  the  "right"  solution  for 
hyperbolic  conservation  laws. 

However,  our  nimwrical  results  indicate  that  (1.1)  has  waak  solutions  with  shock-type 
singularities . 

Acknowledgement!  This  research  was  motivated  by  a  suggestion  of  Professor  John  A.  Nohel. 


2.  ANALYTICAL  THEORY. 

I-  this  Motion,  we  give  •  now  proof  for  the  global  existanca  of  smooth  solutions  in 

the  cut  of  small  data.  (The  local  existence  can  be  shown  using  the  ideas  of  Kato  [9] , 

[10],  cf.  [17]).  The  ideas  used  in  this  proof  will  serve  as  a  guideline  in  our  analysis  of 

the  numerical  scheme  in  section  4.  we  assume  that  the  kernel  b  has  the  form  (1.4)  and 

that  we  are  dealing  with  spatially  periodic  solutions  of  (1.1)>  u(x  +  2*,t)  -  u(x,t).  For 

simplicity,  we  also  aasusa  )  *  1<  (For  4  *  ♦,  a  similar,  but  more  technical  analysis  is 

possible,  sm  the  resuurks  below).  Under  these  conditions,  the  substitution  v  -  ux, 
t  -X.(t  -  s) 

w  “  ut,  gA  -  -ut  ♦  /  e  ♦(ux(x,s))j|de  yields  the  system 


M  M 

(2.1)  w  -  ♦  (v)  -  l  K  g  -  (  l  K  )w  +  f 

x  1-1  i-1 


"  *Vj  "  V  +  J,  V*4  *  W>  '  ' 


This  is  a  perturbation  of  a  hyperbolic  system.  In  order  to  make  it  symmetric  hyperbolic, 

y  _ 

we  definei  w(y)  -  /  /  ♦' (y)  dy,  0(y)  -  »'{«  1 (y ) >  and  set  v  -  w(v).  Then  (2.1) 

0 

assumes  the  form 

v  -  B(v)wx 


(2.2)  w  -  8(v)v  -  l  K  (g  ♦  w)  +  f 

*  i-1  1  1 


M 

“  "V9j  *  *  £  Ki<9i  *  w)  "  f 


Clearly,  0  is  positive  in  a  neighborhood  of  0.  If  +  *  ♦,  (1.1)  can  still  be 

transformed  to  a  syssMtric  hyperbolic  system  after  differentiating  the  equation  [17].  An 
analysis  similar  to  the  following  can  then  be  applied. 
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Th*  analyst*  will  consist  of  thro*  stops ■ 

(i)  Show  that  th*  linear isation  of  (2.2)  at  v  -  0,  w  •  0,  gi  -  0  generates  a 


semigroup  of  negative  type.  As  a  consequence,  the  in homogeneous  initial  value  problem  has 
a  unique  bounded  solution  for  t  e  10,*),  if  th*  inheangeneous  term  is  bounded. 

(ii)  Show  that  th*  same  property  holds  for  the  linearization  of  (2.2)  at 
v  -  vQ(x,t),  w  -  wQ(x,t),  where  vQ ,wQ  are  small  in  an  appropriate  norm. 

(iii)  Os*  a  contraction  argument  for  th*  nonlinear  problem. 

A  ^ 

For  brevity,  let  us  writ*  s  «  (v,w,gf ,. . .gR) .  The  operator  setting  up  th*  linearisation 
of  the  right  hand  side  of  (2.2)  at  '  x  -  0  is  denoted  by  A.  This  operator  acts  Fourier- 

*  r  4  Jrv  *  r  *■  4kx 

componentwise,  i.e.  with  «  “  l  s.  *  we  have  As  ■  L  A.Z.  *  ,  where  the  matrix 

ke*  ka 

A^  is  given  by 


The  characteristic  equation  of  Ak  is 


(2.3) 


l2  +  kV(0) 


.  ,  N  K. 

kV(o>  I 


i-i 


0 


All  eigenvalues  of  A^  have  negative  real  parts,  except  a  double  zero  eigenvalue  for 
k  -  0. 

Proof. 

The  function  defined  by  the  loft  hand  side  of  (2.3)  has  simple  poles  at  X  »  -X^,  and 

therefore  has  zeros  between  these  poles.  This  accounts  for  N  -  1  eigenvalues  of  A^.  A 

further  eigenvalue  lies  between  zero  and  -  min  X  .  For  X  >  0,  th*  left  hand  side  of 

1-1 (1)M  1 

(2.3)  is  always  positive.  If  X  is  complsx,  the  imaginary  part  of  (2.3)  yields 
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(2.4) 


2  Re  A  n  A 


-kV 


(0) 


N 

I  I- 

i-1 


ki  + 


0 


i 

Since  the  eign  of  la  j ia  opposite  to  that  of  la  X,  thla  is  only  possible  if 
Re  1  <  0.  For  k  -  0,  it  is  easy  to  see  that  the  eigenvalues  of  A^  are  sero  (two-fold) 
and  -A^. 

The  presence  of  sero  eigenvalues  is  inconvenient ,  but  we  can  use  a  simple  trick  to  get 

2*  _ 

rid  of  thee,  it  can  be  seen  that  the  integrals  I .  -  /  u  (v)dx  and 
2»  M  1  0 

I,  “  /  (I  t—  (g.  +  w)  -  w)dx  are  invariants  of  the  motion  described  by  (2.2), 

0  i-1  Ai  1 


2v 

provided  that  /  fdx  is  sero.  He  liait  our  attention  to  forces  f  satisfying  this 
0 

condition  and  to  solutions  for  which  the  two  integrals  vanish  initially.  (Physically  this 
aeans  that  the  total  force  and  the  total  mcaaentua  are  sero).  Such  solutions  also  solve  any 
aodified  equation,  in  Which  certain  multiples  of  the  two  integrals  are  added  to  the  right 
hand  side  of  (2.2).  Such  a  modification  replacing  the  double  sero  eigenvalue  by  negative 
eigenvalues  can  be  found.  He  shall  refer  to  the  so  modified  equation  as  (2.2 * )  and  to  the 
modified  linearised  operator  as  A'. 

As  k  ♦  •,  two  eigenvalues  of  A^  have  the  form  iikfS(O)  +  0(1),  the  others 
converge  to  distinct  finite  limits.  Hence  is  diagonalisable  for  large  k.  Let  Tk 
be  a  transformation  matrix  such  that  T~^AyT^  is  diagonal.  The  eigenvectors  of  Ky 
setting  up  the  columns  of  Tk  also  have  limits  as  k  ♦  •,  and  thus  Tk  can  be  normalised 
such  that,  as  k  +  «,  it  converges  to  a  limit  T^.  T^  has  the  form 


As  a  consequence,  both  Tk  and  T^  stay  bounded  as  k  ♦  •».  Therefore,  by  applying  the 
transformation 


kes 


V 


lkx 


Ty  -  I 


kes 


TkV 


ikx 
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-OIU.-  -»T»I »  ts 


f 


the  linearized  operator  In  (2.2*)  assumes  diagonal  form  (up  to,  maybe,  a  finite  dimensional 
perturbation),  and  it  is  clear  that  the  type  of  the  semigroup  generated  by  such  an  operator 
is  determined  by  its  eigenvalues.  As  an  easy  consequence,  we  obtain 
Corollary  2.2. 

For  each  f  e  Hn((0,«)  «  */2  MZtJ**2  )  and  each  *Q  e  Hn(V2W»^H'2),  there  exists  a  unique 

a  A  a  m 

solution  s  of  the  equation  z  -  A’z  *  l,  which  satisfies  the  initial  condition 

A  A 

z(x,0)  “  *q(x). 

Here  Hn  denotes  the  usual  Sobolev  space,  and  n  is  arbitrary.  This  concludes  step 
(i)  of  our  analysis. 

For  abbreviation,  let  us  now  put 

F(v,w,g1,...,g||)  -  (v  -  0(v)wx  +  Ij,w  -  B(v)vx  +  I  *i(«i  +  w)  -  Ij, 

i 

tSj  ♦  ^g^  ♦  l^w  -  l  KA(gA  +  w)  +  I^^vU  -  0),w(t  -  OJ.^tt  -  0))^) 

i 

For  n  >  2,  F  is  a  smooth  mapping  from  Hn([0,«)  x  */2tttJ**2 )  into  | 

H^dO,-)  x  */2**irf*+2)  x  H1*”1  (H/2z*irf**2).  Corollary  2.2  implies  that  the  linearized  3 

i 

mapping  DF(0)  has  an  inverse.  As  a  next  step,  we  show  that  DF  is  invertible  not  only 

J 

at  0,  but  in  a  neighborhood  of  0  (step  (ii))e  1 

/ 

Lemma  2,3. 

if  »°  e  Hn((0,»)  X  (n  >  2)  has  sufficiently  small  norm,  then  DF(z**)  has  an 

inverse  which  is  bounded  as  a  mapping  from  Hn([0,“)  x  n/2n ttl**2)  x  Hn(*/2 z%ilP*2 )  into 
h"< 10, •)  x  M/2z%>tT2u 

Sketch  of  the  Proof i 
2 

He  show  L  -invertlbility.  Estimates  for  the  derivatives  are  then  easily  obtained  by 
successively  differentiating  the  equation.  The  problem  lies  in  the  fact  that  the 
perturbations  resulting  from  the  terms  -0(v)wx  and  -0(v)vx  are  not  bounded  relative 
to  A',  a- id  hence  we  need  a  more  refined  Argument  than  standard  perturbation  theory. 

It  is  well  known  [10]  that,  for  B(v° )  e  H2 (*/2i*i*) ,  the  operator  | 

B  i  (v,w)  ♦  ($(v°)  -  B(0))(wx,vx)  is  a  bounded  perturbation  of  a  skew-adjoint  operator 

B.  Therefore,  the  linearisation  of  (2.2’)  has  the  following  structure 


i 


■7' 


2  2 

where  A*  is  as  above,  B  is  skew-adjoint  and  C  t  L  ♦  L  has  small  norm.  Above,  we 

indicated  how  to  construct  a  transformation  T  *  such  that  A£  is  diagonalized, 

or,  if  it  has  a  degenerate  eigenvalue  for  finite  k,  transformed  to  Jordan  form.  Since 

all  eigenvaluea  are  negative,  we  can  choose  T  such  that  T_1 A'T  is  dissipativet 
*  . 

<s,T  A'Ts)  <  -o(z,s)  for  some  o  >  0.  Moreover,  as  k  ♦  •,  we  have  T,  -  T  ♦  o(r) 

with  Ts  of  the  form  (2.5).  As  a  consequence  T  1BT  -  T~^BT—  is  bounded  (B  is  a  first 

order  differential  operator,  but  *m  -  I  is  of  order  minus  one).  Moreover,  it  is  easy  to 

see  that  T^1 BT—  is  still  skew-adjoint.  It  follows  that  T-1(A*  +  B  +  c)T  is 

dissipative.  Thie  implies  the  lemma. 

Me  have  thus  shown  the  invertibility  of  the  linearisation  in  a  neighborhood  of  0. 

At  this  point,  generalised  implicit  function  theorems  (see  e.g.  [21])  would  be  appli¬ 
cable,  however,  the  quasilinear  nature  of  (2.2*)  admits  a  simpler  argument.  Namely, 

*  A  a  a 

we  can  write  P(s)  in  the  quasilinear  form  p(s)  “  L(s)s,  where 
*  *  *  ~  ^ 

l>(s)z'  «*  DP(0)s*  ♦  (0(0)  -  0(v)  )(w^,v^,0,...,O).  Me  wish  to  solve  the  inhomogeneous 

A  A  a 

problem  F(s)  ■  f.  Since  we  have  juet  shown  that  L(s)  is  invertible  for  small  enough 

*  *  .4  *  A  A 

s,  we  try  the  iteration  *m  -  L  (s<|_1  )f.  if  the  Hn-norma  of  f  and  the  starting 

A  A 

value  z1  are  small,  the  sn  stay  in  a  small  neighborhood  of  zero  in  h".  Moreover, 

we  have 

■  UV‘«M  “>  «V(Vi  -  *«>  ■  *«>>*. 

whence  we  find  the  estimates 


,  -  si  <  Clz  ,  -  z  I  , Is  1 
m+1  m  n-1  m-1  m  n-1  a  n 


If  the  *r®  email  enough,  then  the  sequence  s^  converges  in  Hn_1 ,  and  it  is 

immediate  that  the  limit  oust  be  a  solution.  Me  thus  find 


Theorem  2.4. 

If  f  e  h"([0,-)  k  A/2*Ii ■)  and  v(0),w(0),g,(0),...,g  (0)  e  Hn(V2v«i*) 

*  n 

(n  >  2)  have  sufficiently  small  norms,  then  (2.2')  has  a  solution 

,...,gn)  e  h"([0,-)  x  iy2w*»R**2  which  assumes  the  preecribed 


initial  values 


3.  IAX-WEHDROFF  TYPE  SCHEMES  -  LOCAL  ANALYSIS 

In  this  section  we  do  not  require  the  kernel  b  to  have  the  special  form  (1.4).  It 
is  only  assumed  that  b  is  sufficiently  smooth  and  that  b,b*  are  in  i/  (R+ ) • 

Equation  (1.1)  is  transformed  to  a  system  by  setting 
(3.1)  v  -  ux  ,  w  -  ufc  . 

This  yields 


(3.2) 

with  the  initial  data 


(3.3) 


(b)  wfc  »  ♦(v)^  -  b*+(v)x  +  f(x,t) 

(a)  v(x,0)  -  Ug(x),  x  e  ft 

(b)  w(x,0)  -  u^  (x) ,  x  e  R  . 


The  star  denotes  the  convolution 


(3.4) 


b*p(t) 


t 

/  b(t 


0 


s)p(s)ds  . 


He  discretise  (3.2),  (3.3)  using  a  rectangular  mesh.  The  values  of  the  approximations  to 
v,w  at  the  grid  points  are  denoted  by  a  lower  spatial  and  an  upper  temporal  index 
v"  ,w".  The  following  definition  explains  the  details  of  the  notation. 

Pefinitlon  3.1. 

Let  h  >  0,  k  >  0.  For  n  t  ■  and  i  E  Z  we  set  tn  3  nk,  •  ih.  For 

u  «  (u?)._  —  we  define  the  "spatial"  difference  quotients 

X  lClfncl 


.+  n 


A  u 


n 

i 


2h 


•Q. 


and  the  "temporal"  difference  quotient 


n+1  n 
.♦  n  ui  "  U1 
4  ui - k - * 

The  Idea  of  the  Lax-Wendroff  method  la  to  approximate  v*xj/tn+i )#  w<xi»tn+1 )  by  the 
truncated  Taylor  series 


(3.5) 


(a)  v(Xl,Vl)  -  v(xi.tn)  ♦  kvt(xi<tn)  ♦  r  vtt(xi#tn) 

k2 

(h)  *<VW  -  «<*t.tn>  ♦  ta*t<VV  *  r  wtt<Xl'tn> 


and  to  substitute  for  the  t-derivativea  the  corresponding  express! one  obtained  from 
(3.2).  The  x-derivatives  in  these  expressions  are  approximated  by  spatial  difference 
quotients  in  a  symmetric  way.  This  yields 

2  k 

(a)  v"*1  -  v"  +  kAw"  +  ~  [A+a“*(v")  -  A+A~(b  •  *(v"))  +  Af"] 

(3.6) 

k 

(b)  wf1  -  w"  +  k[A*(v")  -  A(b  •  *(v"))  +  f") 

2  k 

+  —■  (A+(#’(“  (v”+  v"_1))A_w")  -  b(0)A*(v")  -  A(b"  •  «v"))  +  ( f t ) "l 


for  n  >  0,  1  e  Z,  where  we  denoted 

(3.7) 


fl  “  f<VV'  ‘Vi  “  WV 


and  *  is  the  discrete  convolution  operator 

k 


(3.8) 


*  p"  '  k  2  ajnb(tn-j,pj 

J“0  J 


He  require  that  the  weights  a  fulfill 

J« 


(3.9) 


‘  “jn  «  l,  |a},~1  “  V  <  C  * 


where  C  is  independent  of  m.  For  example,  for  the  trapezoidal  rule  cl  «  a  • 

un  nn  t 


o ,  «  1  for  j  -  1  (1)  n  -  1  holds.  Tha  Initial  data  aro 

jn 


(3.10) 


(a)  -  u*^)  ,  1  €  I 


(b)  w”  -  u^)  .  i  «  I 


If  *  Is  at  laast  second  order  accurate  for  sufficiently  smooth  p  and  b,  i.e.  If 

k  . 

(3.11)  |b  *  pn  -  /  b(t  -  s)p( s)ds |  «  0(kZ) 

0 


holds,  then  the  Lax-Hendroff  Method  (3.6).  (3.10)  is  second  order  accurate  (as  far  as  the 
local  discretisation  error  Is  concerned),  provided  that  the  data  and  solution  are  smooth 
enough. 

Thus  our  scheme  is  second-order  consistent.  In  order  to  show  convergence,  we  have  to 
prove  stability  [11].  Since  the  linearisation  is  Lipschits  continuous  with  a  Lipechits 
constant  of  order  o(^) ,  it  is  sufficient  to  show  stability  for  the  linearisation  at  the 
exact  solution,  it  is,  however,  not  enough  to  show  stability  for  the  linearisation  at  sero, 
even  if  we  are  dealing  with  snail  solutions.  Nevertheless,  we  analyse  the  stability  at  the 
trivial  solution  first,  since  this  will  also  give  us  a  guideline  for  dealing  with  the 
variable  coefficient  problem.  He  thus  linearise  (3.6)  at  v™  =  w™  =  0,  which  is  the 
solution  obtained  for  f"  =  0,  uQ (x)  3  u^(x)  =  0.  This  yields  the  constant  coefficient 
scheme 


(a)  “  Ay^  ♦  |  (♦,<0)dVy"i  -  ♦•(O)dV(b  *  y^)]  ♦  g! 


n 

u 


(3.12) 


(b)  «+y"A  -  ♦,<0)Ay"i  -  *M0)A(b  •  y^) 

)( 

♦  |  [♦•(0)A*A"yji  -  b(0)  ♦’  (0 ) Ay™^  -  ♦*«»A(b*  •  y^)]  +  g^ 
for  n  >  0,  i  €  S,  with  the  initial  data 


(3.13) 


(a) 

0 

y1i 

”  y10i  ' 

it! 

(b) 

0 

y2i 

“  y20i  ' 

lex 
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The  stability  analysis  will  proceed  in  the  discrete  L2-space 

‘Vie*1*}  *  "  (h  l  l^lV'2  <  -> 

y 

(Once  stability  in  !>  is  known,  it  is  easy  to  deduce  stability  is  higher  order  discrete 
Sobolev  spaces). 

*•  assume  that  the  initial  data  in  (3.13)  and  the  inhoabogeneous  terms  in  (3.12)  take 
values  in  bjj. 

Since  we  wish  to  apply  Fourier  analysis,  we  extend  (3.12)  to  the  whole  real  line  by 
defining  L  (R)-functions  g™  «  g™  ,  y10  ,  y20  such  that  they  take  the  prescribed  values 

at  the  grid  points.  We  then  regard  (3.12),  (3.13)  as  equations  on  all  of  R.  The  spatial 
difference  operators  are  defined  for  L2(lt) -functions  in  the  obvious  way 


(6*y)(x) 


YiJL-f  h)  -  y(x) 
h 


etc. 


2  * 

For  y  €  I*  (R)  we  denote  the  Fourier  transform  by  yi 

(3.15)  y ( x)  -  /  elx*y(s)ds 

r2t  -m 

Maw  we  prove 
Learna  3.1. 

bet  0  <  |i  -  *  ✓♦•(0)  <  1  and  assuae  that  ♦  ,(0>  >  0,  b,b‘  «  lVr*-),  and  (3.9)  holds. 
Then  the  stability  estimate 


(3.16) 


max 

1-0(1  )n 


2  ♦  2 


)  < 


holds  for  h  sufficiently  small.  Here  y| ,y^ 
Proof,  we  apply  Fourier  transfora  and  set  Yn 
yields  the  equation 


denotes  the  solution  of  (3.12), 


(3.13). 
sh.  This 
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(3.17) 


(3.18) 


y"*1^)  -  B(«)*n(«)  -  *'(o) 


♦MOT  ‘CO*  D\  k  . 


i  •*—  sin  w 
/  ♦MO) 


(b  »  y"(e) ) 


(  °  J(b'  •  y"(a)  +  b(0)y"(a))  ♦  kGn(s) 

✓  ♦'(0)  1 


V  u2  2 

B(«)  *  I  ♦  i  — - A  atn  a  ♦  A  (coi  u  -  1) 

/r"f’7o7  *  51 


(  °  ’)  • 
\  ♦•(0)  0  / 


He  transform  A  to  ita  diagonal  (ora 


(3.19) 


and  aat 


(3.20) 


A  -  1  diag(/V(0 ),-/♦* (0))K_1,  i  .  (  _  _  ) 

\  ✓♦•(0)  -/♦MO)  / 


Vn(a)  -  *-t Tn(B> 


He  thua  obtain  the  new  difference  acha 


k  v"(e)  ♦  v"(a) 


(3.21)  ^’(a)  -  dl.,(,< •),;<«))  -  (SjlJ  I  |)(b  *  (-3 - j-X— )) 


k  v"(B)*w"(B> 


v"(a)«-v"(s) 


-ikK  Lll"  «)(*»•  •  (---_ yj - )  ♦  b(0)  -J—yX - )  J  w"V\s 


(3.22) 


*(«)  ■  1  +  i|i  ain  utp  (coa  a  -  1 ) 


An  eaay  calculation  ahowa  that  |s(a)|  <1,  <u  *  2iw,  ttl,  and  z(2i*)  -  1,  provided 
that  0  <  M  <  1.  In  the  following,  we  only  uae  that  |s(a)|  <  1  ♦  0(k).  Hence,  lc(w))n 
ia  uniformly  bounded  for  n  <  j*.  For  stability  on  hounded  intervals,  we  need  not  be 
concerned  with  0(k) -perturbations  on  the  right  aide  of  (3.21)  (aee  (19),  section  3.9).  It 
ia  thua  sufficient  to  study  the  perturbed  problaa 

(3.23)  h"*’  -  diag(s(a),s(a))Wn  -  $}$}  (f‘“j  I  J)(b  *  p-l-"?.))  -  kH*" 


mm mm* 


which  we  rewrite  as 


k  w"  ♦  w*  n-1 


(3.24)  »"  -  dia,(*"(«) ,;"(«) )W°  -  I  {)(*•  *  (  ♦  *  <»3 


where  the  operator  H  is  defined  as 


H(g3)  -  l  diag(sn"3(ie>,*n"3<*»>G3 

j-0  1-0 


(3.25) 


;?*:3 


Summation  by  parts  yields  (with  u3  «  (~-r — -)) 


(3.26)  l  «n-1'B(»  -  *)(b  *  u")  -  (1  -  *n)(b  *  u"-1) 


"r^  n  » 1 1  n.  .  Vi  .  J  V, 

♦  2  <«  -  s  )  (b  *  u  -b*u) 

s>0 


.  *  "n-1  n-1  "r2  .  *B 

b  u  ••  l1  tkVi,^iV 


♦  *  jlo  <ai.-ib-i-i  '  V-i^1  '  k‘"V.“o 


Moreover,  not.  that  •J#1^b-f1.j  -  «j-b>.J  -  <Bj  -  «jB)  ♦  ■jJVl-i  *  Vj> 

and  that  b|(| ^  j  -  fc^_  j  -  0(k).  Using  (3.9),  one  easily  obtains  froa  (3.26) 


v?  ♦  w"  n-1 


:  >  •  (-4-1))  *<ct„  -x  hi". 

*  '  *-0  j-0(1 )n-1 


with  some  constant  C.  If  we  choose  T  such  that  y  -  CT  <  1,  we  get  from  (3.23)  that 

(3.28)  eax  lw3l  «  y  max  lw3l  +  C.  Iw° I  +  C  l  max  IH3! 

1-0 (1)n  1-0(1 )n-1  l-0(1)n-1 

and  hence 


(3.29) 


■ax  Iw^l  <  .  (c.lw0!  +  c,T  aax  IH^I) 

J*0(1 )n  1-71  i  J*0(1 )n-1 

hold*  for  t  <  T. 
n 

This  argument  can  he  iterated.  He  My  pose  a  new  initial  value  problem  at  the  grid 
point  nearest  to  t  *  T  and  continue  the  solution  into  the  interval  (T,2T)  etc.  This 
yields  an  estiMte  of  the  fora 

(3.30)  aax  Iw^l  <  - - — —  (c  Iw  I  ♦  c.T  mi  IH^I) 

kT<t.<(k+1)T  (1  -  7)  3“0(1 )n 

for  k(l.  (3.16)  iaaediately  follows. 

Ns  see  froa  this  proof  that  solutions  of  the  linearised  problea  can  grow  at  aost 
exponentially.  The  integral  has  been  treated  as  a  ‘lower  order-  perturbation  and  we  Mde 
no  use  of  the  fact  that  it  has  a  daaping  influence,  in  fact,  no  conditions  on  the  sign  of 
the  kernel  were  needed  here.  A  convergence  stateaent  reflecting  the  damping  requires  a 
aore  sophisticated  analysis,  for  a  special  case  this  is  dealt  with  in  section  4. 

As  Mntioned  above,  we  have  to  study  the  linearisation  at  the  exact  solution 
v(x,t) , w( x, t )  of  (3.2),  (3.3)  rather  than  the  linearisation  at  the  trivial  solution.  This 
yields  the  schema 


(3.31) 


(A)  y£’  -  y”A  ♦  kAyjt  ♦  lA*A'(*Mv(xi,tn)y"t)  -  A+  A*(b  •  ♦’(vjy^Jl  ♦  kg" 


1i 


(*>)  y^1  -  y^A  ♦  h[6(V(f)y"l>  -  A(b  •  (♦•(v>y"i>) 


n 

f2i 


♦  (A+(**(“  (v(xt,tB)  ♦  v(*1-1,tn)))A”yJ 

♦  1  ♦'<!  ‘^I'V  +  v,xi-1'tn>,,(yU  *  y1,i-1,4"w<xi'tn,) 


(3.32) 


b(0)A(t'(v)y^1)  -  b*  •  (*,(v)y"i>)  +  kg^ 


(a) 

0 

y1i  _ 

y101 

(b) 

0 

y21  _ 

y201 
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Laws  3.2. 

Assume  0  <  C  <  j*  ux  /V  (»(x,t) )  <  1  hold*  and  assume  that  b,  4,  41  are 
,  +  «e*,te[0,tj 

aaeoth,  b,b'  «  t  (1  I,  f  )  0.  If  the  exact  solution  v,w  is  smooth  on  *  *  (0,  t] ,  than 
an  estimate  of  tha  font  (3.16)  holds  for  tha  solution  of  (3.32),  (3.32)  on  I0,t). 

Proof. 

Again  wo  may  neglect  tarsi  of  ordor  k  on  tha  right  hand  sido  of  (3.31),  which  laavos 
us  with  tha  perturbed  problem 


<•»  <  -  vu +  kiv2i * r  -  f *  <*ivu»  *  *4 


(3.33) 


(b) 


“  4  +  k*Xi  *  r  4V4  -  *“>  :  <4li»  +  kh2 


21 


whore  ♦"  -  ♦*  (vU^tJ ),  *"  -  *,(w<x1,t|l)  ).  With 


(3.34) 


(3.33)  roads 


*2  ♦-  k 

*2 _  a T  a  < 


(3.35) 


'f1  ’  V"  +  kA"4’1  2 


v;  ■  -  v;  ♦  kA*;av?  +  (a?>  VaV"  - 


j-H’lb*  t^)) 


As  before,  wa  diagonal Ire  a"  by  setting 


(3.36) 


h\fi  -m- 


+  khV 


so  that  a"  -  i"j"(b")-1.  wa  substitute 


(3.37) 

thus  obtaining  (after  extending  all  grid  functions  to  real  functions t 


V"  ■  «K 


16- 


rm^  *' 


(3.42) 


Z(x,«)  -  I  ♦  ij"(x)  j  lin  «  ♦  (JB(x))2  ~  ( cox  «  -  1) 
h  h2 


whara  m  has  tha  hm  meaning  aa  bafora.  if  £  /"♦"(*")  «  X  <  1 

n 


for  x  e  K,  than 

ls(x,w)l  <  1  for  x,w  e  X.  It  follows  from  tha  theoreai  on  paga  121  in  [19]  that 


-17- 


(3.43) 


<  1  ♦  Ck 


IL  ' 

"  l  (*)*L  <m) 


•here  C  >  0  can  be  choeen  independently  of  n.k  for  tR  e  (O.T).  The  raat  of  the  proof 
proceeds  in  tha  mb*  lunar  an  for  laaaM  3.1. 

Doing  Kollor'o  (111  nonlinaar  atabliity-conoistanoy  principlo,  wa  obtain  tha  following 
convergence  result t 
Theorem  3.1. 

Lot  b.4,*  bo  ouff iciantly  smooth,  b.b*  «  L1  (•*■),  V  >  0.  Aloo  assume  that  tha 
data  f.Ug.a.  are  aawoth.  Uo  a  eonoaquaneo,  (3.2),  (3.3)  haa  a  aaooth  aolution  v,w). 

/♦' (v(x,t) )  <  1  ho Ida  and 


Chooaa  tha  scheme  (3.0,  (3.10)  auch  that  0  <  e  <  £  wax 

xe»,te[o,  r] 

•uoh  that  (3*9),  (3*11)  ara  aatiafiad#  Than  thara  la  a  uniqua  aolution 


T  ■  *Yi  p  w  •  of  (3.6),  (3*10 )#  which  fulfilla  tha  convarganca 

aatinata 

-  *<VV>i«'  2  +  '<*J  ’  "‘VVW  2  4  ^ 


(3.44) 


whara  D  can  ba  choaan  indapendant  of  h,  tR  <  T. 

*y  •PPlylhD  tha  aaaa  analyaia  to  tha  "differenced"  equation.,  convergence  aatiaataa  in 
higher  Sobolev  nona  can  be  obtained. 

(*  clear  that  tha  nwerical  performance  of  tha  method  vary  much  dapanda  on  how 
"wall*  tha  aecond  order  accurate  integration  rule  (3.8)  integratea  b  and  b’  and  on  how 
(•ff*  darivativaa  of  v  and  w  get.  Particular  care  auat  ba  taken  for  a  claaa  of 
practically  important  kemela  with  tha  following  behavior 


(3.45) 


|bU,(o)|  -  e“0/C,  i  (  I, 


o<e<i,  a  >  o 


Tha  trapacoidal  rule  approximation 

t 

I 

0 


(3.46) 


/  IbMolldo  -  ^  (i  ♦ - 


1  -  a 
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(3.6),  (3.10)  also  requires  thla 


in  unatabla  unlaaa  ~  <  const,  and  tha  axplicit  sch< 
aarioua  restriction. 

In  tha  next  chapter,  we  shall  deal  with  the  special  case  where  the  kernel  is  a  sue  of 
exponentials.  Ha  show  in  this  case  how  stiffness  problems  can  be  avoided,  also,  we 
coapate  the  convolutions  recursively,  which  is  less  tine  consisting  than  using  formula 
(3.8). 


4 


M  imiw  that  tha  karnal  ia  of  tha  fora  (1.4)  and  sat 


t  -X.(t-T) 


(4.1) 


v  -  u 


t**  ■  v*i  “  Kl  1  • 


Than  (1.1)  is  aquivalant  to  tha  systaa 


vt  ’  "x 


(4.2) 


*,  •  ♦<»)_  -  I  «,  ♦  f(x,t> 
x  t*1  * 


•it  ■  V(v,«  -  Vt 


and  tha  initial  conditions  u(x,0)  -  Ug(x),  at(x,0)  -  u.,(x)  bacons 

(4.3)  v(x.O)  “  u*(x),  w(x,0)  -  u^x),  t ^(x.O )  -  0 

Following  tha  racipa  (3.5),  wa  obtain  tha  following  Lax-Wsndroff  discratixation  for  (4.2) 


(4.4)  a^1  -  a"  ♦  kAw*  ♦  uV«(w")  -  A  J  ♦  Af"] 


-r  ■ - 1  «jt  ♦  «"i 


t-i 


♦  r  (4+<*,(i  ♦  Vim'-i>  -  j,  vK> + 1  Vu  + 


H 

I 

1-1 


■u  ■  «U  *  MKtA»(a“)  -  Xt.“4)  ♦ 


♦  §-  (KtA^(»M^  <a"  ♦  vJ_t))A"w")  -  lt(KtA<»(v")  -  X,*^)) 


V1  “  “o^i1'  wl  ”  ui,xi,»  *ti  “  0 


for  1(1,  n  >  0,  subjact  to  tha  initial  conditions 

(4.5) 

Mara  a"  danotas  tha  approximation  to  ▼(x1,t||)  ate.  Tha  last  aquation  of  (4.4)  is  stiff 
for  X(  >  1.  wa  obtain  tha  growth  function 


(4.6) 


*tt(Xtk)  "  1  "  *„*  ♦  \  (Xtk)2  , 
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which  fulfills  |Wj|  <  1  iff  0  <  Xgk  <  2  holds.  Ths  schene  (4.4)  will  thus  bs  unstable 


unless  k  ux  X.  <  2  holds  (for 
i 


rheological  applications  X^  nay  be  as  large  as  10® 


(22]l).  However,  this  can  be  repaired  by  Making  the  schasw  (4.4)  aeei-iaplicit,  e.g.  by 

k^  2  n  V  i 

replacing  the  ten  —  Xji^  by  —  Hjijj  .  In  that  case  we  get 

1  -  X  k 

(4.7)  u  (X  k)  -  - r— * - - 

1  "  2  (X»k> 

and  |u^|  <1  if  0  <  Xjk  <  -1  ♦  /s  or  X^k  >  2.  Moreover,  we  have  0  <  w^  <  1  for 
0  <  Xgk  <1  or  Xgk  >  2.  The  growth  function  of  the  fully  Implicit  scheee  ie  negative 
(but  less  than  one  in  aodulus)  for  X^k  >  1  ♦  /i.  since  this  produces  oscillations  in  the 
ntsaerical  solution,  the  sani-iqplicit  scheee  is  to  be  preferred. 

The  convergence  of  the  schane  (4.4)  is  analysed  in  the  sane  fashion  as  in  section  3, 
the  nontrivial  step  again  being  the  stability  proof.  In  the  present  situation,  stability 
■ust  be  shown  globally  in  tine,  and  not  only  on  finite  tine  Intervals.  We  asaune  ths 
situation  given  in  section  2,  i.e.  ♦  •  solutions  are  spatially  2«-perlodic  (of  course 
the  nesh  else  is  chosen  as  an  integral  divisor  of  the  period,  and  f,uQ,u1  are  snail. 
Rather  than  studying  the  stability  of  (4.4),  we  investigate  the  Lax-Wendroff  discretisation 
of  the  synnetric  hyperbolic  font  (2.2).  Although  the  two  schanes  are  not  equivalent,  they 
only  differ  by  higher  order  terns  negligible  for  the  stability  analysis.  Equation  (2.2) 
leads  to  the  following  schenei 


vi+1  “  v"  *  Wlv")*"  +  f-  [0<v"><A*(0<5  v"  +  i  vt-1,A’vi, 


nt.  n  .  k 


^  a  /A  _ ®  *  A _ ^  v  i" _ ^ * 


A(  I  K^gj^  ♦  w")  +  fj)}  +  ••<»“>•(»“>  <Aa")2] 


*!  "  \  *  k[«(v")4v"  -  l  Kjgj^  ♦  w")  ♦  f"] 


(4.7)  +  —  [6(v")4+(B(i  v"  ♦  -J  v"_, )*-*")  ♦  0,(vJ)*(v")Av"*AwJ 


•  *  V*V«£!l  +  wi)  *  ®(v?,4vi}  +  ftJ 


iH-1  n  ,  ,  r  ,  ,  n  n,  .  p  „  .  n  n.  _n i 

9ii  “  »ti  +  k[-Xt‘»U  +  V  +  J,  V»ai  *  V  -  fJ 


♦  t+Xt(»U  +  "i*  '  +  l  Km{-X„(,!ll  +  wl>  +  Wvi,4wlJ  _  ftJ 


We  first  deal  with  the  stability  of  the  trivial  eolation  v  -  w  -  «  0.  The 

linearisation  of  the  right  hand  aide  at  this  solution  has  the  fora 


Jr  +  -  2 

1  “  I  +  kL.  ♦  kX.,4  +  T~  (L,  A  L.A  +  L  L  A  ♦  L  L  A  ♦  L„) 


"1  1 


0  1  10 


where  Lq,Lj  are  the  constant  coefficient  Matrices 


0(0)  0  ...  0 


0  •  •  •  0 
0 


N 

I  \ 


■X1  *  \  \ 


■XM  +  \  K« 

m-1 


-K- 


-*St 


■A,  ♦  K,  •  •  • 


\ 


+K1  •  •  •  *  *M  / 
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For  the  following,  we  iium 

(4.8)  j“  8(0)  -  const.  <  1  . 

If  K  *  0  is  fixed,  than  L(K)  ■  I  +  k(Lg  +  *3^)  +  0(k2).  Since  the  eigenvalues  of 
Lq  ♦  KL,  have  negative  real  parts  (<-f .  section  2),  the  eigenvalues  of  L(K)  are  inside 
the  unit  circle.  For  K  ■  0,  L(K)  has  a  double  eigenvalue  one,  the  remaining  eigenvalues 
are  inside  the  unit  circle.  The  eigenvalue  one  eay  be  transformed  away  in  a  fashion 
analogous  to  section  2  and  need  not  concern  us  further. 

We  have  to  investigate  the  behaviour  of  L(K)  as  K  ♦  ».  In  this  case,  at  least  one 
of  the  terms  sln(Kh),  cos(Kh)  -  1  is  large  compared  to  h.  The  eigenvalues  of  L(K) 
are  in  first  order  given  by  those  of 

l4l  sin(Kh)L  +  (cos(Kh)  -  1  >L2  . 
n  i  ,4  i 

n 

This  operator  has  an  M-fold  eigenvalue  1  and  the  simple  eigenvalues 

k  it2  2 

Hi  -  8(0  )sin(Kh)  ♦  —  8  (OHcos(Kh)  -  1).  If  (4.8)  holds,  these  two  eigenvalues  lie  on 
h  h2 

an  ellipse  inside  the  unit  circle.  A  simple  perturbation  analysis  shows  that  the  M-fold 
eigenvalue  one  is  perturbed  into  M  distinct  eigenvalues  inside  the  unit  circle,  and  their 
distance  from  the  unit  circle  is  of  order  k.  As  in  section  2,  there  is  a  matrix 
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TOO  such  that  T_1 (K)L(K)T(K)  is  dissipative;  TOC), T-1  (1C)  are  bounded  independently 
of  K,  and  for  K  large  T(K)  has  the  form  (2.5).  Thus  the  trivial  solution  is  stable. 

As  in  section  2,  we  must  show  stability  for  the  variable  coefficient  problem  when  v 
is  in  a  neighborhood  of  0.  There,  we  made  use  of  the  fact  that  the  principal  part  of  the 
differential  operator  was  skew- adjoint.  In  the  discrete  case,  we  have  to  investigate  the 
operator 


L(vi)  -  I  ♦  kL1(vi)A  ♦  vi-t)A" 


where  L^tv)  - 


8(v) 

0 

0 


0  ...  0 

0  ...  0 


All  other  contributions  to  the  linearisation  are  of  order  k.  Let  us  put  A  «  ( v^)  A+. 

Then  the  adjoint  of  A  is  A*  -  -A^L,^)  and  it  is  immediately  verified  that 

k  k2 

L1(vi)  -  I  ♦  -  (A  -  A*>  -  (AA*  +  A*A>  +  O(klvl)  . 

A  simple  calculation  yields 

2 

KI  ♦  |  (A  -  A*)  -  ~  (AA*  ♦  A*A))sl2  -  Isl2  -  i  k2  K  (A  ♦  A*)sl2  +  o(£)3  . 

If  £  is  chosen  sufficiently  small  (but  still  of  order  1),  then  an  argument  similar  to 
section  2  guarantees  the  stability  of  the  variable  coefficient  problem,  we  thus  arrive  at 
the  following  result. 

Theorem  4.1. 

Assume  that  the  assumption  of  section  2  holds  and  that  — 

h 

Then  the  convergence  estimate  (3.44)  holds  uniformly  in  time. 


is  chosen  small  enough. 


5.  NUMERICAL  EXPERIMENTS 

For  all  computations  we  used  a  kernel  b  of  the  fora  (1.4)  and  we  employed  the  scheM 

given  by  (4.4)  and  (4.S)  with  the  semi-implicit  modification.  The  spatial  mesh  sire  was 

prescribed  and  the  temporal  mesh  size  was  determined  at  each  step  such  that  the  stability 

k  / - rT 

conditions  were  satisfied.  In  particular,  we  chose  —  max  /  ♦' ( v. )  <  0.8.  The 

h  i  1 

calculations  were  performed  at  the  VAX  of  the  Mathematics  Research  Center,  University  of 

Ifisconsin-Madison  in  double  precision  arithmetic. 

To  solve  the  initial  value  problem  numerically,  we  introduced  artificial  (far  out) 

boundaries  at  x  *  ±X,  and  the  boundary  conditions  v(iX,t)  -  w(±x,t)  ■  0  (since  the 

initial  data  vanish  at  t**) .  This  introduces  an  additional  error  of  order 

0(  max  | v(±X,t) I  +  max  |w(lX,t> | ).  in  all  the  computations  reported  below,  this 
t«(0,T)  t«(0,T] 

quantity  can  safely  be  neglected.  He  performed  a  convergence  test  for  the  problem  (3.2), 

(3.3)  with  +(v)  -  f(v)  •  v  ♦  {  v3,b(«)  -  0.4e  °,  where  v(x,0),  w(x,0),  f(x,t)  were 

3  2 

2  -r-* 

chosen  such  that  the  problem  has  the  exact  solution  v(x,t)  -  (1  -  x  )e  , 

2 

-f-t 

w( x , t )  -  -xe  .  Table  1  shows  the  errors  ey,  ey  (in  the  discrete  L  -norm)  of  v 
and  w,  rasp,  and  the  corresponding  convergence  rate  given  by  In ("-^y > / in  2  at  two 
different  t-values  and  for  the  maximal  errors  for  t  €  [0,1].  Obviously,  the  scheme  is 
second  order  accurate  tor  this  smooth  solution. 

5 

Table  2  shows  that  the  L  -errors  of  v  and  w  decay  as  t  increases.  The  reason 
for  this  is  the  dissipative  effect  of  the  Volterra  term. 

The  following  calculations  were  done  using 

(5.1)  *(v)  -  *(v)  »  2v  ♦  5v2  +  25v3 

(5.2)  b(0)  -  0.4e"a  +  0.2e-20 
and  the  Initial  data 

(5.3)  (a)  ve(x,0)  -  ev0<x>  -  e(1  -  3x  -  x2  +  x3)e_x  /2 


(b) 


we(x,0) 


ew0(x) 


e(1 


2,  -x2/2 

x  )e 
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with  0  <  e  <  1.  The  force  f (x,t)  wee  eet  equal  to  zero.  4*  Is  strictly  positive  for 
all  vex.  figures  1  and  2  show  Vq(x),  wq(x)  resp. ,  and  the  next  plots  show  t-sectlons 
of  v  and  w,  l.e.  they  show  v  and  w  as  functions  of  x  for  fixed  t- values.  For 
e  “  1  ("large"  Initial  data)  the  dissipative  influence  of  the  Volterra  tens  is  not  strong 
enough  to  avoid  singularities.  Figures  3-5  and  6-8  show  the  evolution  of  v  and  w 
resp.  A  shock-type  singularity  appears  at  t  “  0.057,  x  •  0.4.  Then  a  second 
(t  •  0.073,  x  »  -1.6)  and  a  third  (t  «  0.12,  x  «  2)  shock  develops. 

At  this  point  we  want  to  remind  the  reader  that  the  existence  of  ehocks  for  equation 
(1.1)  has  not  been  proved  yet,  actually  there  is  no  theory  of  weak  solutions  at  all. 
However,  it  has  been  shown  that  smooth  solutions  eay  cease  to  exist  after  a  finite  time 
[8],  because  vx  and  wt  tend  to  infinity. 

Table  3  ehows  the  maximal  values  of  the  difference  quotients  Av",  Aw"  for 
■  1  and  i  <  0  (corresponding  to  the  singularity  at  x  ■  -1.7).  Halving  the  mesh 
size  approximately  doubles  these  values,  which  means  that  the  differences 
v"+1  -  v"  ,  w"+1  -  w"  within  the  shock  are  ptactlcally  independent  of  the  mesh  size. 

Since  the  Lax-Wendroff  method  is  an  artificial  viscosity  method,  we  cannot  expect 
completely  sharp  shocks.  A  shock  layer  of  thickness  0(h)  develops  around  the  shock  ([19], 
section  12.14).  This  is  illustrated  by  Figures  10-14,  which  show  the  left  shock  in  Fig.  9 
for  various  mesh  sizes.  He  have  e  -  ~  and  t  -  0.631.  The  width  of  the  shock  layer  is 
(constantly)  about  3h  (grid  points  are  marked). 

The  "overshoot"  (see  also  Figures  3-8)  is  typical  for  Lax-Wendroff  method  ([19], 
section  12.14)  and  is  due  to  artificial  dispersion.  The  high  wavenumber  components  of  the 
solution  have  a  smaller  wave  speed  and  thus  lag  behind  the  shock  front.  The  width  of  the 
"overshoot  layer"  decreases  with  h.  Outside  the  shock-layer  and  overshoot  region  the 
solutions  coincide  up  to  the  plot  accuracy  for  h  -  0.01,  h  *  0.02,  h  ■  0.04. 

Our  convergence  discussion  in  the  previous  section  does  of  course  not  apply  to 
solutions  with  singularities,  but  it  is  clear  that,  if  the  Lax-Wendroff  method  converges 
boundedly  almost  everywhere,  then  it  converges  to  a  weak  solution  [19]. 


Theref or*  the  praaanted  nuaerical  evidence  indicate*  that  weak  solutions  u  of  (1.1), 
such  that  ufc  and  ux  have  shocks,  exist. 

For  decreasing  e  the  relative  effect  of  dissipation  becoaea  stronger.  For  e  -  ~ 
the  breakdown  of  smooth  solutions  occurs  at  t  »  2.7,  while  the  second  derivatives  of  the 
solution  of  the  corresponding  quasilinear  wave  equation 
<5.4,  utfc  -  ♦(ux,x 

(with  the  saae  initial  data)  blow  up  already  at  t  ■  2.1. 

Figure*  16-21  show  the  evolution  of  v  and  w  for  e  -  ~.  No  singularities  occur 
for  t  «  [0,20],  the  dissipative  eechanisa  of  the  Volterra  ten  seeas  to  produce  globally 
saooth  solutions  her*.  In  the  Figures  1S-21,  L2N  denotss  the  L2-nora  of  v,w  at  the 
given  tin*  t.  The  decay  of  the  L2-noraa  with  increasing  t  is  shown  in  Table  4.  It  is 
clear  froa  Figures  15-21  that  the  L  -noras  also  tend  to  sera  as  t  ♦ 

Figures  21-26  show  the  corresponding  plot  for  the  wave  equation  (5.4)  (without  the 
integral  tera).  The  first  derivatives  of  v  and  w  blow  up  at  t  ■  4.7,  and  the  energy 
given  by  (1.3)  of  course  reaains  constant  for  t  >  0  (apart  froa  artificial  viscosity 
effects  which  show  up  in  the  fourth  digit  of  the  L2-nora). 
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•v(h) 

Rate 

e„(h) 

Rate 

h  -  0.1 

4.7439746 

x  io“3 

2.5703589  x  io'3 

t  •  0.49 

h  -  0.05 

1.2233024 

x  lo“3 

1.96 

5.916561  x  io“4 

2.12 

h  -  0.02S 

3.091415 

x  lo“4 

1.98 

1.432932  x  io-4 

2.05 

h  -  0.1 

6.1822914 

x  10"3 

2.2217676  x  io'3 

t  ■  0.96 

h  -  0.05 

1.5S63728 

x  10" 3 

2.00 

5.859864  x  io'4 

1.92 

h  -  0.025 

3.883662 

x  io-4 

2.00 

1.553665  x  io'4 

1.92 

h  -  0.1 

6.198451 

x  10"3 

3.3348472  x  io'3 

MX 

h  -  0.05 

1.5611899 

x  10-3 

1.99 

7.615021  x  io-4 

2.13 

teio.ij 

h  -  0.025 

3.906411 

x  IQ-4 

2.00 

1.809902  x  io-4 

2.07 

Table  1.  Errors  and  Convergence  Rates 
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Table 


t 

s^Oi  «  0.1) 

ew(h  -  0 

.2) 

1. 

6.0477392  a 

10“3 

2.6280316 

a  10"3 

3. 

3.3230997  a 

10"3 

2.6443695 

a  10-3 

5. 

2.1187529  a 

10"3 

1 .9864532 

a  10-3 

7. 

1.4750333  a 

10"3 

1.4025019 

a  10“3 

9. 

1.0429657  a 

10'3 

9.730021 

a  10“4 5 

11. 

6.102893  a 

10-4 

5.120333 

a  10“4 

13. 

2.335227  a 

10-4 

1.883220 

a  10“4 

15. 

2.232786  a 

10-4 

1.348227 

a  10"4 

Table  2. 

Decay  of 

Errors 

Av 

da 

h  ■  0.1 

8.748 

53.19 

h  •  0.05 

17.90 

113.8 

h  -  0.025 

34.41 

223.6 

3.  Numerically  obtained  value*  for  vx,wx  at  t  •  1,  x  ■ 


t  L2N  :  V 

0.  5.3925  *  10-2 

5.  1.20131  *  10-2 

10.  3.7581  *  10"3 

15.  1.5276  *  10"3 


L2N  i  N 
2.88243  *  10"2 
1.60351  »  10"2 
4.81279  a  10‘3 
1.8037  *  10“3 


Table  4.  Decay  of  L,2.norm#  due  to  dissipation 
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0.00000 


Figure  3.  Development  of  shocks  for  large  data. 


Development  of  Shocks  for  large  data. 
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Figure  5.  Development  of  shocks  for  large  data 
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Figure  6.  Development  of  shocks  for  large  data 
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Figure  7.  Development  of  shocks  for  large  data 


'iqure  8.  Development  of  shocks  for  large  data 
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Figure  9.  V(X,t  =  0.62099)  for 


Figure  11.  Dependence  of  shock  width  on 


Figure  13.  Dependence  of  shock  width  on 


6.46861  L2N=  0.0085771 


Figure  18.  'Small'  Data, 
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.92462 


5.59847  L2N=  0.04084971 
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.73299  L2N=  0.05871420 


5.59847  L2N=  0.05866761 
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